﻿using System;
using System.Collections.Generic;
using System.Text;
using System.ComponentModel;


namespace KorTransLib
{
    public class Jeodezi
    {
        const double deg2rad = Math.PI / 180;
        const double rad2deg = 180.0 / Math.PI;

        /// <summary>
        /// Estimates the distance between corFrom and corTo.
        /// </summary>
        /// <param name="corFrom">Coordinate parameter.</param>
        /// <param name="corTo">Coordinate parameter.</param>
        /// <returns></returns>
        public static double distVincenty(Coordinate corFrom, Coordinate corTo)
        {
            TypeConverter corcon = TypeDescriptor.GetConverter(typeof(DecimalDegreeCoordinate));
            DecimalDegreeCoordinate ddCoord1 = (DecimalDegreeCoordinate)corcon.ConvertFrom(corFrom);
            DecimalDegreeCoordinate ddCoord2 = (DecimalDegreeCoordinate)corcon.ConvertFrom(corTo);
            
            double a = 6378137, b = 6356752.3142, f = 1 / 298.257223563;  // WGS-84 ellipsiod
            double L = (ddCoord1.lon - ddCoord2.lon) * deg2rad;
            double U1 = Math.Atan((1 - f) * Math.Tan(ddCoord1.lat * deg2rad));
            double U2 = Math.Atan((1 - f) * Math.Tan(ddCoord2.lat * deg2rad));
            double sinU1 = Math.Sin(U1), cosU1 = Math.Cos(U1);
            double sinU2 = Math.Sin(U2), cosU2 = Math.Cos(U2);

            double lambda = L, lambdaP, iterLimit = 20;
            double cosSqAlpha, cos2SigmaM;
            double sinSigma;
            double cosSigma;
            double sigma;
            do
            {
                double sinLambda = Math.Sin(lambda), cosLambda = Math.Cos(lambda);
                sinSigma = Math.Sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) +
                  (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));
                if (sinSigma == 0) return 0;  // co-incident points
                cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
                sigma = Math.Atan2(sinSigma, cosSigma);
                var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
                cosSqAlpha = 1 - sinAlpha * sinAlpha;
                cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
                if (Math.Abs(cosSqAlpha) < 1e-12)
                    cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)
                var C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
                lambdaP = lambda;
                lambda = L + (1 - C) * f * sinAlpha *
                  (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
            } while (Math.Abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0);

            if (iterLimit == 0) return 0;  // formula failed to converge

            var uSq = cosSqAlpha * (a * a - b * b) / (b * b);
            var A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
            var B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
            var deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
              B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
            var s = b * A * (sigma - deltaSigma);

            return s;
        }

        public static void googleTest()
        {
            double lat = 38.441523;
            double lon = 27.110984;
            int zoom = 12;

            double radLat = lat/180*Math.PI;
            double radLon = lon/180*Math.PI;

            double mercx = radLon;
            double mercy = Math.Log(Math.Tan(radLat)+1/Math.Cos(radLat));

            double cartx = mercx+Math.PI;
            double carty= Math.PI-mercy;
            double pixelxz0 = cartx *256/(2*Math.PI);
            double pixelyz0 = carty *256/(2*Math.PI);
            int tilex = (int)(pixelxz0 *Math.Pow(2,zoom)/256);
            int tiley = (int)(pixelyz0 *Math.Pow(2,zoom)/256);
            int serverZoom = 17 - zoom;
        }
    }
}
